// -----------------------------------------------------------------------------
//
//                            procedure dsinit
//
//   this procedure provides deep space contributions to mean motion dot due
//     to geopotential resonance with half day and one day orbits.
//
// Author:
//   Jeff Beck
//   beckja@alumni.lehigh.edu
//   1.0 (aug 7, 2006) - update for paper dav
// original comments from Vallado C++ version:
//   author        : david vallado                  719-573-2600   28 jun 2005
//
//   inputs        :
//     Cosim, Sinim-
//     Emsq        - Eccentricity squared
//     Argpo       - Argument of Perigee
//     S1, S2, S3, S4, S5      -
//     Ss1, Ss2, Ss3, Ss4, Ss5 -
//     Sz1, Sz3, Sz11, Sz13, Sz21, Sz23, Sz31, Sz33 -
//     T           - Time
//     Tc          -
//     GSTo        - Greenwich sidereal time                   rad
//     Mo          - Mean Anomaly
//     MDot        - Mean Anomaly dot (rate)
//     No          - Mean Motion
//     nodeo       - right ascension of ascending node
//     nodeDot     - right ascension of ascending node dot (rate)
//     XPIDOT      -
//     Z1, Z3, Z11, Z13, Z21, Z23, Z31, Z33 -
//     Eccm        - Eccentricity
//     Argpm       - Argument of perigee
//     Inclm       - Inclination
//     Mm          - Mean Anomaly
//     Xn          - Mean Motion
//     nodem       - right ascension of ascending node
//
//   outputs       :
//     em          - eccentricity
//     argpm       - argument of perigee
//     inclm       - inclination
//     mm          - mean anomaly
//     nm          - mean motion
//     nodem       - right ascension of ascending node
//     irez        - flag for resonance           0-none, 1-one day, 2-half day
//     atime       -
//     d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433    -
//     dedt        -
//     didt        -
//     dmdt        -
//     dndt        -
//     dnodt       -
//     domdt       -
//     del1, del2, del3        -
//     Ses  , Sghl , Sghs , Sgs  , Shl  , Shs  , Sis  , Sls
//     theta       -
//     xfact       -
//     xlamo       -
//     xli         -
//     xni
//
//   locals        :
//     ainv2       -
//     aonv        -
//     cosisq      -
//     eoc         -
//     f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543  -
//     g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533  -
//     sini2       -
//     temp        -
//     temp1       -
//     theta       -
//     xno2        -
//
//   coupling      :
//     getgravconst
//
//   references    :
//     hoots, roehrich, norad spacetrack report #3 1980
//     hoots, norad spacetrack report #6 1986
//     hoots, schumacher and glover 2004
//     vallado, crawford, hujsak, kelso  2006
//  ----------------------------------------------------------------------------*/

function dsinit(
    cosim,  emsq,   argpo,  s1,     s2,     s3,     s4,
    s5,     sinim,  ss1,    ss2,    ss3,    ss4,    ss5,
    sz1,    sz3,    sz11,   sz13,   sz21,   sz23,   sz31,
    sz33,   t,      tc,     gsto,   mo,     mdot,   no,
    nodeo, nodedot,       xpidot, z1,     z3,     z11,
    z13,    z21,    z23,    z31,    z33,    em,     argpm,
    inclm,  mm,     nm,     nodem, ecco,   eccsq) {

    // /* --------------------- local variables ------------------------ */
    var twopi = 2.0 * Math.PI,
    aonv  = 0.0,
    q22    = 1.7891679e-6,
    q31    = 2.1460748e-6,
    q33    = 2.2123015e-7,
    root22 = 1.7891679e-6,
    root44 = 7.3636953e-9,
    root54 = 2.1765803e-9,
    rptim  = 4.37526908801129966e-3,
    root32 = 3.7393792e-7,
    root52 = 1.1428639e-7,
    x2o3   = 2.0 / 3.0,
    znl    = 1.5835218e-4,
    zns    = 1.19459e-5,
    // sgp4fix identify constants and allow alternate values
    irez = 0,
    // TODO: FIGURE OUT HOW TO GET GLOBALS
    //global tumin mu radiusearthkm xke j2 j3 j4 j3oj2
    //tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2,
    // /* -------------------- deep space initialization ------------ */
    d2201 = 0,
    d2211 = 0,
    d3210 = 0,
    d3222 = 0,
    d4410 = 0,
    d4422 = 0,
    d5220 = 0,
    d5232 = 0,
    d5421 = 0,
    d5433 = 0,
    del1  = 0,
    del2  = 0,
    del3  = 0,
    atime = 0,
    xfact = 0,
    xlamo = 0,
    xli   = 0,
    xni   = 0,
    // variable defined later
    ses, sis, sls, sghs, shs, sgs, sghl,
    dedt, didt, dmdt, sdhl, shll, domdt, dnodt, dndt, theta,
    cosisq, emo, emsqo, eoc,
    g201, g211, g310, g322, g410, g422, g520, g533, g521, g532,
    sini2, f220, f221, f321, f322, f441, f442, f522, f523, f542, f543,
    xno2, ainv2, temp1, temp,
    g200, g300, f311, f330;

    if ((nm < 0.0052359877) && (nm > 0.0034906585)) {
        irez = 1;
    }
    if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) {
        irez = 2;
    }

    // /* ------------------------ do solar terms ------------------- */
    ses  =  ss1 * zns * ss5;
    sis  =  ss2 * zns * (sz11 + sz13);
    sls  = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
    sghs =  ss4 * zns * (sz31 + sz33 - 6.0);
    shs  = -zns * ss2 * (sz21 + sz23);
    //   // sgp4fix for 180 deg incl
    if ((inclm < 5.2359877e-2) || (inclm > Math.PI - 5.2359877e-2)) {
        shs = 0.0;
    }
    if (sinim !== 0.0) {         // WARN unreliable floating point comparison
        shs = shs / sinim;
    }
    sgs  = sghs - cosim * shs;

    // /* ------------------------- do lunar terms ------------------ */
    dedt = ses + s1 * znl * s5;
    didt = sis + s2 * znl * (z11 + z13);
    dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
    sghl = s4 * znl * (z31 + z33 - 6.0);
    shll = -znl * s2 * (z21 + z23);
    //   // sgp4fix for 180 deg incl
    if ((inclm < 5.2359877e-2) || (inclm > Math.PI - 5.2359877e-2)) {
        shll = 0.0;
    }
    domdt = sgs + sghl;
    dnodt = shs;
    if (sinim !== 0.0) {        // WARN unrealiable
        domdt = domdt - cosim / sinim * shll;
        dnodt = dnodt + shll / sinim;
    }

    // /* ----------- calculate deep space resonance effects -------- */
    dndt   = 0.0;
    theta  = (gsto + tc * rptim) % twopi;
    em     = em + dedt * t;
    inclm  = inclm + didt * t;
    argpm  = argpm + domdt * t;
    nodem  = nodem + dnodt * t;
    mm     = mm + dmdt * t;
    // //   sgp4fix for negative inclinations
    // //   the following if statement should be commented out
    // //if (inclm < 0.0)
    // //  {
    // //    inclm  = -inclm;
    // //    argpm  = argpm - pi;
    // //    nodem = nodem + pi;
    // //  }

    //  /* - update resonances : numerical (euler-maclaurin) integration - */
    //  /* ------------------------- epoch restart ----------------------  */
    //  //   sgp4fix for propagator problems
    //  //   the following integration works for negative time steps and periods
    //  //   the specific changes are unknown because the original code was so convoluted

    // /* -------------- initialize the resonance terms ------------- */
    if (irez !== 0) {
        aonv = Math.pow(nm / xke, x2o3); // TODO GLOBAL xke

        // /* ---------- geopotential resonance for 12 hour orbits ------ */
        if (irez === 2) {
            cosisq = cosim * cosim;
            emo    = em;
            em     = ecco;
            emsqo  = emsq;
            emsq   = eccsq;
            eoc    = em * emsq;
            g201   = -0.306 - (em - 0.64) * 0.440;

            if (em <= 0.65) {
                g211 =    3.616  -  13.2470 * em +  16.2900 * emsq;
                g310 =  -19.302  + 117.3900 * em - 228.4190 * emsq +  156.5910 * eoc;
                g322 =  -18.9068 + 109.7927 * em - 214.6334 * emsq +  146.5816 * eoc;
                g410 =  -41.122  + 242.6940 * em - 471.0940 * emsq +  313.9530 * eoc;
                g422 = -146.407  + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
                g520 = -532.114  + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;
            }
            else {
                g211 =   -72.099 +   331.819 * em -   508.738 * emsq +   266.724 * eoc;
                g310 =  -346.844 +  1582.851 * em -  2415.925 * emsq +  1246.113 * eoc;
                g322 =  -342.585 +  1554.908 * em -  2366.899 * emsq +  1215.972 * eoc;
                g410 = -1052.797 +  4758.686 * em -  7193.992 * emsq +  3651.957 * eoc;
                g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
                if (em > 0.715) {
                    g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
                }
                else {
                    g520 =  1464.74 -  4664.75 * em +  3763.64 * emsq;
                }
            }
            if (em < 0.7) {
                g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21  * eoc;
                g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
                g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4  * eoc;
            }
            else {
                g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
                g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
                g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;
            }
////
            sini2 =  sinim * sinim;
            f220 =  0.75 * (1.0 + 2.0 * cosim + cosisq);
            f221 =  1.5 * sini2;
            f321 =  1.875 * sinim  *  (1.0 - 2.0 * cosim - 3.0 * cosisq);
            f322 = -1.875 * sinim  *  (1.0 + 2.0 * cosim - 3.0 * cosisq);
            f441 = 35.0 * sini2 * f220;
            f442 = 39.3750 * sini2 * sini2;
            f522 =  9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) +
                                       0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq));
            f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim +
                                                  10.0 * cosisq) + 6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq));
            f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim + cosisq *
                                       (-12.0 + 8.0 * cosim + 10.0 * cosisq));
            f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim + cosisq *
                                       (12.0 + 8.0 * cosim - 10.0 * cosisq));
            xno2  =  nm * nm;
            ainv2 =  aonv * aonv;
            temp1 =  3.0 * xno2 * ainv2;
            temp  =  temp1 * root22;
            d2201 =  temp * f220 * g201;
            d2211 =  temp * f221 * g211;
            temp1 =  temp1 * aonv;
            temp  =  temp1 * root32;
            d3210 =  temp * f321 * g310;
            d3222 =  temp * f322 * g322;
            temp1 =  temp1 * aonv;
            temp  =  2.0 * temp1 * root44;
            d4410 =  temp * f441 * g410;
            d4422 =  temp * f442 * g422;
            temp1 =  temp1 * aonv;
            temp  =  temp1 * root52;
            d5220 =  temp * f522 * g520;
            d5232 =  temp * f523 * g532;
            temp  =  2.0 * temp1 * root54;
            d5421 =  temp * f542 * g521;
            d5433 =  temp * f543 * g533;
            xlamo =  (mo + nodeo + nodeo - theta - theta) % twopi;
            xfact =  mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - no;
            em    = emo;
            emsq  = emsqo;
        }

        // /* ---------------- synchronous resonance terms -------------- */
        if (irez === 1) {
            g200  = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
            g310  = 1.0 + 2.0 * emsq;
            g300  = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
            f220  = 0.75 * (1.0 + cosim) * (1.0 + cosim);
            f311  = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
            f330  = 1.0 + cosim;
            f330  = 1.875 * f330 * f330 * f330;
            del1  = 3.0 * nm * nm * aonv * aonv;
            del2  = 2.0 * del1 * f220 * g200 * q22;
            del3  = 3.0 * del1 * f330 * g300 * q33 * aonv;
            del1  = del1 * f311 * g310 * q31 * aonv;
            xlamo = (mo + nodeo + argpo - theta) % twopi;
            xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no;
        }

        // /* ------------ for sgp4, initialize the integrator ---------- */
        xli   = xlamo;
        xni   = no;
        atime = 0.0;
        nm    = no + dndt;
    }

    return [em,     argpm,  inclm,  mm,     nm,     nodem, irez,
            atime,  d2201,  d2211,  d3210,  d3222,  d4410,  d4422,
            d5220,  d5232,  d5421,  d5433,  dedt,   didt,   dmdt,
            dndt,   dnodt,  domdt,  del1,   del2,   del3,   xfact,
            xlamo,  xli,    xni];
}